#!/bin/bash

# averaging across time points to correct for Fuzzy Ripples
3dTstat -mean -prefix INV1_AP.nii.gz -overwrite moco_Basis_TI1_100.nii
3dTstat -mean -prefix INV1_PA.nii.gz -overwrite moco_Basis_TI1_101.nii

3dTstat -mean -prefix INV2_AP.nii.gz -overwrite moco_Basis_TI2_100.nii
3dTstat -mean -prefix INV2_PA.nii.gz -overwrite moco_Basis_TI2_101.nii

# Taking the ratio of inversion times to corret for Rx bias fields,
# proton density, and T2*
3dcalc                                                                       \
    -a          INV1_AP.nii.gz                                               \
    -b          INV2_AP.nii.gz                                               \
    -expr       'min(a/b,1)'                                                 \
    -overwrite                                                               \
    -prefix     UNI_AP.nii
3dcalc                                                                       \
    -a          INV1_PA.nii.gz                                               \
    -b          INV2_PA.nii.gz                                               \
    -expr       'min(a/b,1)'                                                 \
    -overwrite                                                               \
    -prefix     UNI_PA.nii

# removal of noise outside the brain
LN_MP2RAGE_DNOISE                                                            \
    -INV1    INV1_AP.nii.gz                                                  \
    -INV2    INV2_AP.nii.gz                                                  \
    -UNI     UNI_AP.nii.gz                                                   \
    -beta    20                                                              \
    -output  UNI_AP_denoised.nii
LN_MP2RAGE_DNOISE                                                            \
    -INV1    INV1_PA.nii.gz                                                  \
    -INV2    INV2_PA.nii.gz                                                  \
    -UNI     UNI_PA.nii.gz                                                   \
    -beta    20                                                              \
    -output  UNI_PA_denoised.nii

# Bias field correction in SPM
/Applications/MATLAB_R2024a.app/bin/matlab                                   \
    -nodesktop                                                               \
    -nosplash                                                                \
    -r          "Bias_field_script_job"

# Distortion estimation
# [PT] **** add -noXdis and -noZdis to have AP-PA ONLY
# so, probably have to _check orientation_ OR resample before/after
3dQwarp                                                                      \
    -plusminus                                                               \
    -pmNAMES    Rev For                                                      \
    -pblur      0.05 0.05                                                    \
    -blur        -1 -1                                                       \
    -noweight                                                                \
    -minpatch   9                                                            \
    -source     mUNI_PA_denoised.nii                                         \
    -base       mUNI_AP_denoised.nii                                         \
    -overwrite                                                               \
    -prefix     blip_warp.nii
        
# let's look at the distortion field in the phase encoding direction
# just for QA
# [PT] this step might not be necessary IF we restrict the warp axes above
# for the phase encode dist corr
3dcalc                                                                       \
    -a          blip_warp_For_WARP.nii'[1]'                                  \
    -prefix     warpfield_phaseencode.nii                                    \
    -expr       'a'                                                          \
    -overwrite

# help to reduce fat ghosts --- maybe replace this with MIN(...) 
3dcalc                                                                       \
    -overwrite                                                               \
    -prefix     unwarpedmean.nii                                             \
    -a blip_warp_For.nii                                                     \
    -b blip_warp_Rev.nii                                                     \
    -expr 'mean(a,b)'                                                         

# Echo spacing of T1234 is 1.29 ms, segmentation and inplane GRAPPA is
# 14.  This makes the effective echos pacing (relevant for
# distortions) 0.09214 ms

# Echo spacing of functional protocol is 1.12 ms, segmentation and
# inplane GRAPPA is 4.  This makes the effective echo spacing 0.28 ms.

# Estimation of the distortion field of functional data
# [PT] put some of these numbers UP ABOVE bc the user will NEED to precalc 
# these
3dcalc                                                                       \
    -a          blip_warp_For_WARP.nii                                       \
    -expr       'a*((-1)*0.28/0.09214)'                                      \
    -prefix     warptofunctional.nii                                         \
    -overwrite

# Applying the distortion field to the ANATOMICAL so it will sit nicely with
# the EPIs
# [PT] 
3dNwarpApply                                                                 \
    -source     unwarpedmean.nii                                             \
    -master     unwarpedmean.nii                                             \
    -prefix     mached2functional.nii                                        \
    -nwarp      'warptofunctional.nii'                                       \
    -overwrite

